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ABSTRACT 

The Taiwanese-American Occultation Survey (TAOS) aims to detect serendipitous occultations of 
stars by small (~1 km diameter) objects in the Kuiper Belt and beyond. Such events are very rare 
(< 10~^ events per star per year) and short in duration (~200 ms), so many stars must be monitored 
at a high readout cadence. TAOS monitors typically ^500 stars simultaneously at a 5 Hz readout 
cadence with four telescopes located at Lulin Observatory in central Taiwan. In this paper, we report 
the results of the search for small Kuiper Belt Objects (KBOs) in seven years of data. No occultation 
events were found, resulting in a 95% c.l. upper limit on the slope of the faint end of the KBO size 
distribution oi q — 3.34 to 3.82, depending on the surface density at the break in the size distribution 
at a diameter of about 90 km. 

Subject headings: Occultations, Kuiper belt: general, Comets: general, Planets and satellites: forma- 
tion 



1. INTRODUCTION 



The size distribution of Kuiper 
been accurately measured down 
p > 30 km (Fuentc s ct al. 2009; 
2008t iFraser fc KavelaarsI I2008t 



Belt Objects has 
to diameters of 
Fuentcs fc Holmiffli 
iFraser et all 12005 



Bernstein et al.l 12004 iLuu fc JewittI I2002D. while 



Eraser fc KavelaarT i 200S ) have provided measurements 
down to diameters l) > 15 km. The size distribution of 
large KBOs approximately follows a power law of the 
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down to a diameter of I? « 90 km. A clear break 
in the size distribution has been detected near 
90 km (iFuentes et al.l 120091: IFraser fc KavelaarsI 120091: 
iBernstein et al.l I2004D . giving way to a shallower dis- 
tribution for smaller objects. Various models of the 
formation of the Kuiper Belt have been developed 
dBenavide z fc Carnpo Bagati n 2009: Kenvon & Bromlej 
2009; Pan & Sari|_j2005l Kenvon & Bromley 2004 
20011: i Bcnz & Asphaug 1999; Kcnyon & Luu 1999a. 
Davis fc Farinella,1997; ■Stern.J996i : .Duncan et al...l99£ 



which predict this break in the size distribution, but 
these models make vastly different predictions on the 
size distribution for objects smaller than the break 
diameter. These models are based on different scenarios 
of the dynamical evolution of the Solar System and the 
internal structure of the KBOs themselves. It is gen- 
erally assumed that the outward migration of Neptune 
stirred up the orbits of objects in the Kuiper Belt, which 
subsequently underwent a process of coUisional erosion, 
giving rise to a shallower size distribution. The simplest 
of these models also predict a small-end power law size 
distribution of the form 



dn 
dD 



cx D- 



(2) 



where the slope q depends primarily on the internal 
strength of the KBOs. Solid objects held together by 
material strength give rise to larger slopes, while the 
weaker gravitationally bound rubble piles are more eas- 
ily broken up and give rise to smaller values of q. A 
measurement of the size distribution of smaller objects 
(down to D > 100 m) is needed in order to constrain 
these models. Furthermore, such a measurement would 
provide import ant infor mation on the origin of short- 
period comets (iVolk fc Malhotr al 120081: lTa,ncredi et al.l 
2006; Duncan & Levison 1997; Levison fc DuncanI 119*971 
Morbidcni.1997; .Holman fc Wisdom.l993i) . 
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The direct detection of such objects is difficuh be- 
cause they are extremely faint, with typical mag- 
nitudes R > 28, and are thus invisible to sur- 
veys using even the largest telescopes. However, 
a small KBO will induce a detectable drop in the 
brightness of a distant star when it pa sses across 

the line of sig ht to the star |Chang ct alj 120121: 

Schlichting et ali 120121: [Wang et all 12010/ .Bianco et al l 



20101: iSchlichting et alJ 120091: iBickertonet alJ 1 200 



Wang et all 12009^ 'Bianco et al.ri2009t iLiu et al? '2OO8' 



Zhang et all j2008: Bickcrto iTet all 12008^ Ichang et al 
2007; Nihci ct al. 2007; Chang ct al. 2006; "R ogues et aTl 



2006, .2003; ,Coorav 2003; 
Rogues fc Moncuaueti I2OOC 



Coorav &: Farmen 
jBrown fc Websten 



2003 



1997 



Rogues et all 119871: lBailevlll976[ ). High speed photom 



etry is needed to detect these events, given that we ex- 
pect a typical event duration of about 200 ms. Detection 
of such events is further complicated by the fact that 
the sizes of the objects of interest are on the order of 
the Fresnel scale (about 1.4 km for our median observed 
wavelength of about 600 nm and a typical object distance 
of 43 AU), and the resulting lightcurves exhibit signifi- 
cant diffration features (Nihci ct al. 2007). The goal of 
the TAOS project is to detect such occultation events 
and measure the size distribution of Kuiper Belt Objects 
(KBOs) with diameters 0.5 km < D < 30 km. 

To date, 18 events consistent with occultations by 
objects in the outer Solar Syst em have been detect ed. 
The three detections claimed bv I Rogues eiTall l|2006^. if 
they are occultations, are unlikely to have been caused 
by KBOs (the distances are inconsi stent). On the 
other hand, two detections reported bv ISchlichting et all 
(|2009D and ISchlichting et al.l ()2012D are consistent with 
occultation events by 1 km diameter KBOs at 45 AU. 
[Schlichting ct al. (2012) found one additional event con- 
sistent with a KBO occultation, but they reject this event 
because of its high inclination (it was detected at an 
ecliptic latit ude of b = 81.5°). The remaining 12 pos- 
sible events (jChang et al.l 120071) were found in archival 
RXTE measurements of Sco X-1. However, the interpre- 
tation of th ese detections as occultation e vents has been 
challenged (|Jones et al.ll2008l : IBlocker et al..,2009l 

TAOS has been in operation s ince February 2005. The 
system is described in detail in iLehner et al.l ()2009D . In 
summary, the survey operates four 50 cm telescopes at 
Lulin Observatory in central Taiwan. The telescopes are 
very fast (F/1.9), with a field of view of 3 deg^. Each 
telescope is equipped with a camera which utilizes a sin- 
gle 2kx2k CCD imager. All four telescopes monitor the 
same stars simultaneously at a readout cadence of 5 Hz. 
We require that any events be detected simultaneously 
in each of th e telescopes in orde r to minimize the false 
positive rate (jLehner et al.l l2010'). 

High speed cadence is achieved using a special CCD 
readout mode called zipper mode, which is described in 
detail in ILehner et al.l ()2009() . To summarize, instead of 
taking repeated exposures, we leave the shutter open for 
the duration of a data run. At a cadence of 5 Hz, we 
read out a subimage comprising 76 rows, which we have 
termed a rowblock. When the rowblock is read out, the 
remaining photoelectrons on the CCD are shifted down 
by 76 rows as well. If one considers starting with a row- 
block at the far edge of the CCD, photoelectrons are 
collected and then shifted to the next rowblock. Photo- 



electrons are then collected in addition to those from the 
first rowblock, and then they are shifted to the next. By 
the time the original photoelectons are read out, photo- 
electrons have been collected from every star on the focal 
plane in a single rowblock, although at different epochs. 

While zipper mode readout allows us to image every 
star in the focal plane, there are two characteristics of 
the data that reduce the signal-to- noise ratio (SNR). The 
first is that sky background is collected as a rowblock is 
shifted across the focal plane, while photoelectrons from 
a single star are only collected at one location. Further- 
more, it takes about 1.3 ms to read out a single row, 
corresponding to about 95 ms for an entire rowblock. So 
the actual exposure of a star is only 105 ms, while sky 
background is collected for a total of 5.4 s. Second, flux 
is also collected from the stars during the readout stage, 
and the brighter stars will leave significant numbers of 
photoelectrons in the pixels as they are being shifted. 
This flux will appear as bright vertical streaks in the 
images (see the upper left panel of Figure[T]). Neverthe- 
less, zipper mode is successful in collecting high SNR 
photometric data on enough stars for the project to suc- 
cessfully probe a significant fraction of the allowed size 
distribution parameter space. 

In this paper we report on the analysis of seven years 
of zipper mode data collected by the TAOS system. In 
^ we describe the data set used in this analysis. In ^JSj 
we provide a detailed discussion of the analysis pipeline. 
In 21 present the results of the analysis of the data 
set, and finally in Sj5l we discuss our future plans for the 
survey. 

Throughout the remainder of this paper, we use the fol- 
lowing definitions. A data run is a series of high cadence 
multi-telescope measurements on a single field, typically 
for 1.5 hours at a time. A lightcurve is a time series 
of photometric measurements of a single star from one 
telescope in one data run, and a lightcurve set is a set 
of multi-telescope photometric measurements of a sin- 
gle star during one data run. A flux tuple is a set of 
multi-telescope measurements (with a minimum of three 
telescopes) of one star at a single epoch. A lightcurve set 
is thus a time series of flux tuples. 

2. THE SEVEN YEAR DATASET 

Earlier results f rom the TAOS sear ch fo r small KBOs 
were r eported in iZhang et al.l (|2008f l and iBianco et al.l 
(jMol ). These two papers describe results obtained 
from operati n g thr ee telescopes. The dataset used in 
IBianco et all (|2010D ended on 2008 August 2, when the 
fourth telescope became operational. The addition of the 
fourth telescope to the array was beneficial in two ways. 
First, using four telescopes increases the significance of 
any candidate occultation ev ents, wh ile the false positive 
rate remains constant (Lehner et al.l[2010.) . Second, we 
have experienced frequent problems with the reliability 
of our cameras and telescopes. With four telescopes, we 
can still operate with three telescopes when one of the 
systems is shut down for repair. Given that a minimum 
of three telescop es is necessary to ke ep the false positive 
rate low enough (jLehner et alllMoh . we can still collect 
useful data in the event that there is a problem with one 
of the systems. The overall operational efficiency of the 
survey thus improved significantly since we could still 
collect data when one telescope was offline. 
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Table 1 

Dataset parameters. 



Zhang et al. (2008 'I Bianco et al. (2010) This Work 

Start Date 
End Date 
Data Runs 
Light-curve Sets'' 
Total Exposure (star— hours)'' 
Photometric Measurements'' 
Flux Tuples" 



2005 Feb 7 2005 Feb 7 2005 Feb 7 

2006 Dec 31 2008 Aug 2 2011 Sep 8 

156 414 858 

110,554 366,083 835,732 (209,130) 

152,787 500,339 1,159,651 (292,514) 

7.8 X 10^ 2.7 X 10^0 6.7 x 10^° (2.1 x 10^°) 

2.6 x 10^ 9.0 X lO'-* 2.1 X 10^° (5.3 x 10^) 



Values for four-telescope data shown in parentheses. 



I 



Figure 1. Top left: Horizontal subsections of two consecutive rowblocks with a bright moving object in the center. (The boundary 
between the rowblocks is visible as a dark horizontal line bisecting the image. This is a slight excess in dark current from the edge of 
the CCD.) The vertical streaks (see text) from the brighter stars are also evident. Top right: The same rowblock subsections after the 
background streak subtraction. While the streaks from the bright stars are removed, the moving object itself is not completely removed, 
and the background flux near the moving object is over-subtracted. Bottom left: grayscale image representation of the background flux 
subtracted from the rowblocks near those shown in the top panels. Here each row represents the flux subtracted from a single rowblock, 
while the columns correspond to those in the top panels. The vertical streaks in this image are the streaks left by the stars during the 
zipper mode readout. The moving object in original data cause the background to be over-subtracted, and the over-subtracted flux clearly 
shows the trajectory of the moving object across the focal plane. Bottom right: The same background image from the bottom left panel 
after applying the high pass mean and variance filters to remove the streaks from the bright stars. The over-subtracted background flux 
from the moving object is clearly evident. 



The current dataset is sum marized in T able [H 
with the datasets used in iZhang et al.l 



along 

^ ^ and 

iBianco et al.l ()2010[ ). The current dataset is mo r e tha n 
2.3 times larger than that used in I Bianco et al.l (|2010[ ). 

90% of the photometric data in this set was acquired 
within 6° of the ecliptic in order to maximize the event 
rate. The results quoted in this paper are thus applica- 
ble to both the hot (inclination z > 5°) and cold {i < 5°) 
KBO populations. 

3. DATA ANALYSIS 

The data analysis took place in three stages: photo- 
metric reduction of raw image data, the search for oc- 
cultation events, and the detection efficiency simulation. 
T hese stages roughly fo llow the same procedures outlined 
in IBianco et al.l ( 20101 ) , but several improvements were 
made to the pipeline in this round of data analysis. The 
new pipeline is described in the following subsections. 

3.1. Photometric Reduction 

For the first step in the pipelin e, a custom photometric 
reduction pipeline ()Zhang et al.l l20091 is used to measure 
the brightness of each star in the field at each epoch, and 
the resulting measurements are assembled into a time se- 
ries lightcurve for each star. For every star, the lightcurve 
from each telescope is combined into a lightcurve set. 
Next, bad data points (e.g. measurements where the 
star lies near the edge o f the focal plane) are removed 
from the lightcurves (see IZhang et al.ll2009l for details). 

A significant improvement added to this stage in the 
process is the automated flagging of photometric mea- 



surements affected by bright objects (such as satellites 
or airpla nes) moving across the field of view. As dis- 
cussed in IBianco et al.l ()2010t ). several candidate events 
were found in the previous dataset. Visual inspection 
of the relevant images revealed that these events were 
caused by brigh moving objects passing near the star 
that appeared to be occulted. The new algorithm to de- 
tect these bright moving objects is described in i jS.l.ll 

3.1.1. Bright Moving Object Detection 

As mentioned in i )3.11 bright objects moving across the 
field of view during zipper mode image collection can give 
rise to a large number of false positive events. This is 
caused by the background subtraction in the photometric 
reduction. As discussed in ^ the zipper mode readout 
has the disadvantage that the brighter stars leave streaks 
along the columns in the images due to the finite time it 
takes to read out a row from the CCD. These streaks are 
removed during the photometric reduction by subtract- 
ing an estimate of the mode of the column from each pixel 
in the column, wher e the mode is calcu lated for each zip- 
per mode rowblock (jZhang et al.ll2009|) . A bright moving 
object passing over a particular column causes the mode 
of that column to be over-estimated, and the background 
is thus over-subtracted, causing an artificial drop in the 
measured brightness of any star whose aperture includes 
that column, (see the top panels of Figured]). 

In order to detect these events, we save the 3o--clipped 
column mean for each rowblock and column on all four 
cameras. These background data are stored in FITS for- 
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Figure 2. Histograms of p- values for the Pearson's test (left panel) and hypergeometric test (right panel) for each of the data runs in 
the data set. 



mat, and an example is shown in the bottom left panel 
of Figure[TJ Note that each column in this image cor- 
responds to an actual column on the CCD, but the row 
corresponds to a rowblock in the zipper mode data. The 
pixel value is simply the column mean (after 3cr clipping) 
for that rowblock and column. The streaks from the 
bright stars are evident in the image. Also clearly visible 
is the object shown in the top panels moving across the 
image. We note that what appears as a moving object 
is actually the high value for the column average due to 
the moving object. Nevertheless, it is straightforward to 
identify which columns in which rowblocks are affected 
by the over-estimated background. 

First, we normalize the image by applying the same 
rolling clipped mean and variance filter to each of the 
columns we use to remove the slowly varying trends along 
the lightcurves (jLehner et al.ll2010r ). That is, we replace 
each column average value 6mri, where m is the row block 
index and n is the column, with a normalized value dmn , 
defined as 



dmn, — 



(3) 



where fimn is the local (3cr-clipped) mean within column 
n and amn the local variance within column n, calculated 
with window sizes of 21 rowblocks. Note that cr„m is 
calculated after the mean has been subtracted from every 
rowblock background value along the column. 

The bottom right panel of Figure[T] shows the resulting 
values of dmn after application of this filter to the im- 
age shown in the bot tom left panel. We then apply an 
algorithm ()Lutdll980l )^ to find objects in the normalized 
background image. We define an object as either a set 
of 20 or more adjoining pixels with values of dmn larger 
than 1.5, or as a set of 5 or more pixels with values larger 

^ This is the s ame algorithm used in SExtractor 
IIBertin He Arnoutslll99d) to detect objects. 



than 2.5. We then flag any photometric measurements 
where a column affected by a moving object lies within 
the aperture used in the photometry. The removal of all 
flagged data points results in the elimination of all false 
positive events caused by moving objects. 

3.2. Event Search 
3.2.1. Rank Product Statistic 

Our event selection algorithm uses the rank product 
statistic, which is described in detail in iLehner et al.l 
(201_0). This statistic takes advantage of the multi- 
telescope data to give an exact value for the signiflcance 
of any candidate event, even if the underlying distri- 
bution of the flux measurements is not known. Given 
that TAOS samples at a 5 Hz readout cadence, we can- 
not resolve the features of any candidate events in our 
lightcurves. Any occultation event we detect would typ- 
ically appear as a small drop in the flux for one or two 
consecutive measurements (although occultation events 
by larger objects could contain up to five or six consecu- 
tively measured flux drops). We therefore need to have 
an accurate estimate of our false positive rate in order to 
report a credible result. 

To search for events in the data, we calculate the rank 
product statistic on each tuple in the entire data set. To 
calculate this statistic, we take a time series of fluxes 
/i, /2, . . . , /atp (where A'p is the number of points on the 
time series), and replace each measurement /-,- with its 
rank rj. The lowest measurement in the time series will 
be given a rank of 1, and the brightest measurement will 
be given a rank of A^p. For a given lightcurve set, we thus 
replace the time series of flux tuples with a time series of 
rank tuples of the form 

{rij,...,rTj), 

where T is the number of telescopes used. We do not 
know the underlying distribution of the flux tuples in 
each lightcurve set, but, assuming the lightcurves meet 



TAOS Seven Year Results 



5 



—\ T" 



1 km 




_l L_ 



-20 



20 



time (s) 




-5 '- 



H h 



H h 



H h 



3 km 

^ h 



H h 



H h 



H h 



H h 



H h 



-20 



20 



time (s) 



Figure 3. Lightcurves with simulated occultation events from 1 km (top) and 3 km (bottom) diameter objects added in, before and after 
application of the moving average. In each panel, the top lightcurve is the simulated event, the second lightcurve has no moving average 
applied, and the next three have moving averages with window sizes 2, 3, and 5 respectively. 

rank product would be significantly smaller than average, 
leading to a large value of our test statistic z. We set a 
threshold zt for event detection based on choosing an 
upper limit on the expected false positive rate, such that 



the correct requirements (discussed in H3.2.2p . we do 
know the exact distribution of rank tuples, and we can 
thus calculate the significance of any candidate event ex- 
actly. 

The rank product statistic is calculated as 



InTT^. 



P{Z > Zt) = 



(4) 



0.25 



(5) 



If one assumes that the rank distribution for a single 
lightcurve is continuous, then this statistic follows the 
gamma distribution. Given that the ranks are discrete, 
the approximation fails for large values of z, however, the 
probability distribution can be calculated exactly using 
simple combinatorics. 

In the case of an occultation event, one expects the flux 
to drop below the nominal value simultaneously in each 
of the telescopes. The corresponding rank tuple would 
have lower values for all of the ranks, and the resulting 



where Z is a random variable chosen from the rank prod- 
uct probability distribution, 0.25 is our limit on the ex- 
pected number of false positive events in the data set, 
and -/Vt is the total number of tuples used in the event 
search. 

3.2.2. Lightcurve Filtering and Data Quality Checks 

The rank product test statistic follows the expected 
distribution for an individual lightcurve set if and only if 
the following conditions are satisfied (Cochlo 2010). 

• The distribution of measurements in each 



6 



Zhang et al. 



_ 1 1 1 1 1 1 

■1 

' : 

- 1 


1 1 1 1 1 1 _ 

f 

8 km J 


- — ' 'l ' ' ' 

^ 'z — 1 1 1 1 1 


1 1 1 1 

1 1 1 1 


hJ'--/--AvJ\A^ 

5 ^ 1 




1 1 1 1 

5 ^ , , 1 


1 1 1 1 


1 1 1 1 1 1 1 1 1 1 _ 

, 1 , , , 1 , , , 1 , ," 


-20 ( 

time 

_ 1 1 1 1 1 1 


) 20 

1 1 1 1 1 1 _ 


1 : 1 

- 1 


J 30 km : 


: — 1 1 1 1 1 

^ : — 1 1 1 1 1 


1 1 1 j\ 1 1-^ 

1 1 1 1 1— 


^ : — 1 1 1 1 1 




5 - 

1 1 1 1 1 1 


1 1 1 1 


^ - , , 1 , , , 1 , , , 1 



-20 

time (s) 

Figure 4. Same as Figure|3l but for 8 km (top) and 30 km (bottom) diameter objects. 
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lightcurve in the lightcurve set is stationary. 

• Each Hghtcurve in the hghtcurve set is ergodic in 
mean. 

• Each hghtcurve in the hghtcurve set is independent 
of the other hghtcurves in the hghtcurve set. 

However, none of the hghtcurve sets in our data set meet 
these requirements. Changes in the transparency of the 
atmosphere throughout the course of our typical 1.5 hour 
data runs induce fluctuations in the hghtcurves, render- 
ing them non-stationary. Furthermore, these fluctua- 
tions are correlated between the different telescopes, so 
the hghtcurves are not independent of each other. We 
therefore apply a rolling mean and variance filter to each 
point in the hghtcurves, thus replacing each flux mea- 
surement fj with hj, given by 



^3 



(6) 



where /ij is the local mean and Uj the local variance at 
point J, calculated with window sizes of 33 and 151 re- 
spectively. Note that the variance is calculated after the 
rolling mean was subtracted, as was done in Equation|31 
In most cases, this filter removes all of the slowly vary- 
ing trends, and the resulting filtered hghtcurves h meet 
the requirements for the application of the rank product 
statistic. 

However, there are still some data runs where fast 
moving cirrus clouds can induce simultaneous variations 
in the hghtcurves that are not removed by the high 
pass filter. We hav e therefore devised a pair of tests 
(jLehner et al.l 120100 to identify lightcurve sets where 
the filtered hghtcurves are not independent (and likely 
not stationary as well, given that the fluctuations in 
the hghtcurves from the cirrus clouds are not constant 
over time). The tests are based on the fact that the 
rank tuples should be distributed uniformly over the T- 
dimensional rank space (recall that T is the number of 
telescopes used). For the first test, we divide the T- 
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Figure 5. Top panel: histogram of SNR values for lightcurve sets 
in the entire data set. (For each lightcurve set, the minimum SNR 
among all telescopes is used for the histogram.) Lowrer panels: SNR 
values of recovered events from our efficiency simulation for 1 km, 
8 km, and 30 km objects (solid histogram). Also shown, using the 
right vertical axis, is the detection efficiency vs SNR (dotted lines). 

dimensional rank space uniformly into a T-dimensional 
grid. If the lightcurves are independent, we expect that 
each cell in the grid would contain the same number of 
rank tuples. We thus calculate the Pearson's statistic 
on the number of tuples in each cell of the grid. For the 
second test, we only look at the number of rank tuples 
in the lowest ranked cell in each dimension. When look- 
ing at all of the lightcurve sets in a given data run, the 
results from the first test should follow the distribu- 
tion, and the results from the second test should follow 
the hypergeometric distribution. (However, this is only 
true if the lightcurves are independent and identically 
distributed, which is a stricter requirement than station- 
ary. We have thus devised a hlockwise bootstrap method 
to calculate the expected distributions. We still call the 
tests the Pearson's test and the hypergeometric test.) 
We can thus calculate a p-value for both tests for each 
data run to see how well the two test statistics match the 
expected distributions. Histograms of these p- values are 
shown in Figure[5]for each test. If the statistics match the 
distribution, we would expect a uniform distribution of 
p- values in the range of [0, 1]. The large spikes evident 
in the lowest bins of each histogram are the data runs 
that have correlated fluctuations in the lightcurves that 
are not removed by the high pass filter. This leads to 
a non-uniform distribution of rank tuples, meaning that 
the test statistics for each of the lightcurve sets does not 
follow the expected distribution. We thus remove every 
data set with a p-value less than 0.01. Note that this 
means we would only expect to remove 1% of our data 
runs due to random chance. 

In order to perform these tests, we require enough high 
(SNR stars to accurately calculate these two test statis- 
tics. We thus set a minimum of 10 stars with SNR > 10 
in the data run. Any data run that does not meet this 
requirement is excluded from further analysis. Visual in- 
spection of lightcurves in such data runs indicate that 
there are very few stars, and the SNR is compromised, 
usually by very bright sky during a full moon. These 
data runs are not expected to contribute significantly to 
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Figure 6. Plots of fioff vs D for the optimum combinations of 
moving average window size s and SNR cu t s. Th e values are scaled 
by SIbiOi which is Qgff from lBianco et at] 1 12010 ). in order to high- 
light the differences. 

our effective sky coverage, so we lose little by excluding 
them from the final results. Furthermore, as will be dis- 
cussed in !j3.2.4[ most of the stars in these data runs will 
end up being excluded anyway. 

3.2.3. Detection of Multipoint Occultation Events 

The rank product statistic described in §3.2.11 only 
works on one rank tuple at a time. This has the disad- 
vantage of being inefficient at detection of longer dura- 
tion occultation events (such as objects with D > 5 km, 
scattered disk or Oort Cloud objects beyond 100 AU, 
or events measured at large o p positi on angles). How- 
ever, as was shown in iCoehlol (|2010l ) and discussed in 
iLehner et al.l (|2010l ). if we replace a stationary lightcurve 
with a moving average of the form 



1 n 

Gj = —[rii + 
w 



''j+w — l}^ 



(7) 



where hj is the lightcurve after high pass filtering and 
w is the rolling average window size, then the resulting 
lightcurve will also be stationary, and the rank product 
statistic will still be applicable. By applying a moving 
average, we can increase the significance of any candidate 
event since the average lightcurve will exhibit a more 
significant drop in the measured flux when the averaging 
window is centered on the event. 

This is illustrated in Figures |3] and 21 which shows the 
results of the application of a moving average on actual 
lightcurves with simulated events added in. We have 
used three different window sizes w = 2, 3, and 5, in 
order to optimize our detection efficiency. We also looked 
at the lightcurves with no average applied, which we will 
refer to as using a window size w = 1 in order to simplify 
the discussion. For the 1 km diameter simulated event, 
we detect the event using w = 1, 2, and 3. For such 
a small object, using w = 5 actually smoothes out the 
signal, while amplifying other small fluctuations, so that 
the event does not pass the cuts. For the 3 km and 
8 km diameter objects, the events are detected only with 
w = 1 and 2. For the 30 km object, which is longer in 
duration than the other simulated events, the signal is 
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Figure 7. Q^ff vs D using our final selection criteria of ui = 1 
and 3 and SNR > 3. A lso sliown for comp arison are the pre- 
vious TAOS result s from [B lanco et al.l l l20TOtl (dotted line) and 
IZhane et all ll2008l) (dashed line). 

not detectable for w = 1, but it is detected with w — 2, 
3, and 5. 

It is clear that the different window sizes increase the 
sensitivity to different size objects. However, while the 
effective sky coverage of the survey can be increased by 
using a wide array of window sizes, this w ill also increase 
the false positive rate ()Lehner et al.ll2010l ) . Running each 
multiple window size on each lightcurve set is not an in- 
dependent test, given that events can be detected using 
more than one window size, as demonstrated in FiguresO 
andlH It follows that any false positive event could be de- 
tected multiple times using different window sizes. How- 
ever, given that we do not know the underlying distribu- 
tion of flux values in the data set, we can not estimate 
the rate at which this may occur. Therefore, instead of 
estimating the true false positive rate, we can only set an 
upper limit on the false positive rate by assuming that 
searching for events in a lightcurve set after applying the 
moving average is a completely independent test for each 
value oi w. If we have a total of window sizes, we 
thus modify Equation[5] as 



P{Z > zt) 



0.25 

NtN^ 



(8) 



So the more window sizes that are used, the tighter the 
detection threshold must be. The optimization of the 
selection of the window sizes used in this analysis will be 

discussed in 21 

Note that the data quality tests described in §3.2.21 
must be reapplied f or each window size w. As shown in 
iLehner et al.l ()2010t ). application of the moving average 
can highlight small correlations between the telescopes 
that are insignificant in the unaveraged data. There- 
fore using larger window sizes will require removing more 
data runs from the analysis. In cases where we use more 
than one window size, we always choose to be conser- 
vative and remove the union of data sets that would be 
rejected using each value of w. 



The final cut we apply to our data set is to set a thresh- 
old on SNR of each lightcurve set. The goal of this cut 
is to exclude stars which provide little sensitivity to oc- 
cultation events, yet still contribute to the total number 
of tuples Nt, which tightens the selection threshold (see 
EquationH]) . This is illustrated in Figure[5l which shows 
the total distribution of SNR values for every lightcurve 
set in the data set, as well as the SNR distributions of 
stars for which simulated events are recovered for 1 km, 
8 km, and 30 km objects. It is clear that there is little 
contribution to our effective sky coverage for SNR < 5. 
Also shown is the detection efhciency for the same di- 
ameter objects vs. SNR. As expected, the detection ef- 
ficiency drops off at low SNR values, especially for the 
smaller objects. 

The exact value of SNR to use as a threshold depends 
on many factors, especially the set of moving average 
window sizes that will be used. This optimization is dis- 
cussed in 21 

4. RESULTS 

We searched through our data set, using moving aver- 
age window sizes of 1, 2, 3, and 5, (and every combination 
thereof) as well as SNR cuts of 0, 1, 2, 3, 4, 5, 6, and 7. 
For each combination of window size and SNR thresh- 
old, the detection threshold zt was calculated as shown 
in Equation[8l In all cases, no events were found in the 
data set. So the next step in the analysis was to measure 
our effective sky coverage and optimize the selection of 
moving average window sizes and SNR cut. 

The total number of events expected by the survey if 
given by 

dn 



dD 



n,s{D) dD, 



(9) 



where dn/dD is the differential size distribution (the 
number of objects per deg~^ per km along the ecliptic), 
and fieff is the effective sky coverage of the survey. We 
estimate ilcff by implanting a simulated event into each 
lightcurve set and seeing if it is recovered by our selection 
criteria. For this simulation, we choose a set of diameters 
between 0.5 km and 30 km, and for each diameter, we 
implant exactly one event into each lightcurve set. We 
assume every object is at a geocentric distance of 43 AU, 
since the lightcurve shape does not depend critically on 
the distance for most objects (40 to 46 AU) in the Kuiper 
Belt. The tests for each diameter are performed indepen- 
dently, so there is never more than one event implanted 
into any lightcurve set at any time. The effective sky 
coverage is calculated as 



EH w,oi „ 



(10) 



where the sum is over only those lightcurve sets I where 
the added event was recovered, Ei is the length of the 
lightcurve set in time, v^qi is the relative transverse veloc- 
ity between the Earth and the KBO, A is the geocentric 
distance to the KBO (again, fixed to a constant value of 
43 AU), and H is the event cross section. We estimate 
H as 



3.2.4. Signal to Noise Ratio Cuts 



H{D, A, 9,) w f(2V3i^i + oi 



+ 0.A, 



(11) 
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Table 2 

Results from using different window sizes and SNR cuts. 



Window Sizes 


SNR Cut 


Tuples 


2t 
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n^s{D = 1km) 




noff(i? = 30 km) 




1 & 2 


3 


9.49 X 10^ 


1.32 X IQ- 


11 


3.82 


(3.590 ±0.053) X 10" 


7 


(2.788 ± 0.011) X 10- 


5 


1 & 2 


4 


7.13 X 10'' 


1.75 X IQ- 


11 


3.81 


(3.703 ± 0.054) X 10' 


7 


(2.577 ±0.010) X 10- 


5 


1 & 2 


5 


5.63 X 10^ 


2.22 X IQ- 


11 


3.82 


(3.806 ±0.055) X 10' 


7 


(2.277 ±0.010) X 10- 


5 


1 & 3 


3 


9.14 X lO'' 


1.37 X IQ- 


11 


3.82 


(3.467 ± 0.053) x 10" 


7 


(5.065 ± 0.014) X 10- 


5 


1 & 3 


4 


6.87 X lO'' 


1.82 X IQ- 


11 


3.82 


(3.559 ±0.053) x 10" 


7 


(4.430 ± 0.013) X 10- 


5 


1 & 3 


5 


5.42 X lO'' 


2.31 X IQ- 


11 


3.82 


(3.657 ±0.054) x 10" 


7 


(3.770 ± 0.012) X 10- 


5 


1 & 5 


3 


8.60 X 10^ 


1.45 X IQ- 


11 


3.84 


(3.213 ±0.051) X 10- 


7 


(5.594 ± 0.015) X 10- 


5 


1 & 5 


4 


6.46 X lO'' 


1.94 X IQ- 


11 


3.84 


(3.301 ±0.052) X 10" 


7 


(4.824 ± 0.014) X 10- 


5 


1 & 5 


5 


5.08 X lO'^ 


2.46 X IQ- 


11 


3.84 


(3.413 ±0.052) X 10" 


7 


(4.071 ± 0.013) X 10- 
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Figure 8. Solid line: 95% c.l. upper limit on cumulative surface density vs. diameter D (bottom axis) and R magnitude (top axis, 
assuming an alb edo of 4% and a distance of 43 AU) from current TAOS data set. Dotted lines: 95% c.l. upper an d lower limits on surface 
density i^rom Schlichting et al.l f2012). Cross and error bar (S12): Surface density reported by ISchlichtin g et al.l f2012). (Note that this 
is not strictly model independent because it is based on a power law model. However , the result docs not depend very strongly on the 
assumed slo pe.) Solid squares (RXTE): upper limit (right point) reported bv lLiu et al.l 12008) and the best fit surface density (left point) 
(2007) from an occultation search through RXTE observation of Sco X-1. 



reported by [Chang_et_ al . _ , . 

reported by I Wang et al . (2010) using Pan-STARRS guider images. Empty triangle (BKW ): u pper limit reported bv lBickerton et al.l I I200SI ) 



Solid trian gles (PS): upper limits 



using the 1.8 m telescope at DAO. Empty circles (MMT): upper limits reported by 'Bianco et al. (2009) from analysis of trailed images 
obtained with the MMT. Upside down empty triangle (R06): upper limit reported by Roqucs ct al. (2006) from high speed imaging using 
the 4.2 m Herschel Telescope. Empty square ( HST): u pper limit reported by 'Bernstein et al. (2004) from a direct survey using the HST. 
Sol id circle (FGH): surface de nsity reported bv lFuentes et al. (2009) from a direct search using Subaru. Star (FK): surface density reported 
by 'Eraser fc Ka velaarsI I I2009I) from a direct survey using Subaru. (Note that Eras er fc KavelaarsI 120091) assume an albedo of 6% and a 
distance of 35 AU, so they claim that the point plotted at R = 27.875 co rrespon ds to a diameter D = 15 km.) L ines and box labeled 
CB, P, and SD: require d surface density for the Classical Belt ILevison fc Duncanlll997l 'l. Plutinos lMorbidellilll997l ). and Scattered Disk 
IVolk fc Malhotrall200l) respectively to be the source for the observed distribution of Jupiter family comets. 
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where 0* is the angular size of the target star and F is 
the Fresnel scale, which is included to account for the 
minimum cross s ection of an event due to diffraction 
(|Nihei et al.l[2nn7h . 

To optimize the selection of window sizes and SNR 
cut, we looked at the resulting values of fioff for every 
possible combination. There is no clear optimum set of 
cuts, so we decided to optimize our selection based on 
three criteria. 

• Minimize the upper limit of the slope q at the small 
end of the size distribution, assuming a power law 
size distribution anchored at the break diameter 
of 90 km and a cumulative su rface density at the 
break diameter of 5.4 deg" ^ (jFraser fc Kavelaari 
[20091 : TSchlichting et"aI1 120091) . 

• Maximize the se ns itivity at D = 1 km, where 
iSchlichting et al.l (|2009l) and iSchlichting et al.l 
(120121 )" claim the detection of two KBOs. 

• Maximize the sensitivity at 13 = 30 km, in order to 
bring our upper limit as close as possible to the cur- 
rent direct detection limits (iBernstein et al]l2004l : 
iFuentes et aLll2009l : HVaser fc Kavelaarsll2009fr 

The best results came from using cuts on SNR of 3, 4, 
and 5. The results are summarized in Table [2] We found 
that the upper limit on the slope q did not vary signif- 
icantly as long as we limited ourselves to two window 
sizes, and included w = 1. The combination of w = 1 
and 2 and a minimum of SNR > 3 gave the largest effec- 
tive sky coverage at 1 km, however, this was particularly 
bad at 30 km. The best result at 30 km came from w — 1 
and 5 with SNR > 3, but this combination was not very 
good for 1 km. The combinations that worked best for 
both diameters was using w — I and 3 with SNR cuts of 
3 and 4. The resulting values of flcs for these four com- 
binations are shown in FigurelHl In the end, we opted for 
using w = 1 and 3, with a threshold of SNR > 3 as this 
is better at 30 km, and not much worse at 1 km. 

Given our final set of cuts, the resulting plot of flcs vs 
D is shown Figure[7l The points on the curve indicate the 
diameters where we calculated QcS using our detection 
efficiency simulation, and the curve itself is a cubic spline 
fit to these values. 

To set a model independent upper limit, with no de- 
tected events we can eliminate any model which would 
lead us to expect more that three detected events at the 
95% c.l. Our model independent upper limit is shown in 
Figure[51 along with results from other occultation sur- 
veys and direct searches. We note that our upper limit 
at J> = 1 km i s a factor of 4.5 below that reported by 
ISchlichting et a l. (2012), but it is also a factor of 6.7 
higher than their lower limit based upon their two re- 
ported events. Therefore even though we found no events 
at 1 km (despite the fact that our efficiency simulation 
for 1 km diameter objects showed that we are in fact 
sensitive to events similar to those detected by the HST 
survey), our limit is thus still consistent with their pub- 
lished result. 

We also calculate upper limits on the small KBO size 
distribution assuming it follows a power law anchored at 
the detected break at 13 = 90 km. We use Equation[9] 
to calculate the number of expected events as a function 




Figure 9. Number of expected events vs slope q. Solid line: num- 
ber of expected events using a surface density of A'^o^go km = 
5.4 deg~^, as reported by ll^'raser &: Kavela ars' f2009D. Dashed 
line: number of expected events vs. slope q, using a surface den- 
sity o f Af£)>gokm = 38.0 deg~^, as reported bv iFuentes et aLl 
120091 ). Since no events were found, any model predicting more 
than 3 events (dotted line) is excluded at the 95% c.l. 

of slope q. The results are shown in Figure^ We use 
two values for the surface density at D = 90 km, which 
normalize the size distribu tion for smaller diam e ters. Us- 
ing the value reported by iFraser fc Kavelaari (|2009( ) of 
Nd>90 km = 5.4 deg~^ y ields a 95% c.l. upper limit of 
q — 3.82. In comparison, ISchlichting et al.l (|2009[ ) report 
a slope oi q ~ 3.8 ±0.2 based on their claim of two detec- 
tions at D = 1 km, assuming the inclination distribution 
of small KBOs follows tha t of the larger (D > 100 km) 
objects (|Elliot et al.ll2005[ ). The second v alue we use is 
Np^gp km = 38.0 deg~^ as reported by IFuentes et al.l 
(|2009l ). This resuhs in a 95% c.l. upper limit on q of 
3.34. 

5. FUTURE PLANS 

After seven years of observations, it would only be of 
marginal value to continue the survey in its present form. 
So at the end of the data set discussed in this paper, we 
shut down the system for a camera upgrade. Each new 
camera, manufactured by Spectral Instruments, utilizes 
a CCD47-20 frame transfer CCD from e2v. With the 
new cameras, we can image with a readout cadence of 
just under 10 Hz if we use 2x2 binning. The CCDs 
are Ikxlk with 13 /im pixels, as opposed to the 2kx2k 
imagers with 13.5 /im pixels used to collect the current 
data set, so 77% of the field of view is lost. However, 
by moving away from zipper mode to full frame imaging, 
our SNR increases significantly and we estimate that our 
limiting magnitude will increase from R — 13.5 to 15. 
Wc thus expect to be able to monitor a similar number 
of stars. Given the higher readout cadence, we become 
significantly more sensitive to objects with D < 1 km. 
Observations with these new cameras began in November 
2012. 
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